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Abstract. The cometary meteoroid ejection model of Jones and Brown (1996b) was 
used to simulate ejection from comets 55P/Tempel-Tuttle during the last 12 
revolutions, and the last 9 apparitions of 109P/Swift-Tuttle. Using cometary 
ephemerides generated by the Jet Propulsion Laboratory’s (JPL) HORIZONS Solar 
System Data and Ephemeris Computation Service, two independent ejection schemes 
were simulated. In the first case, ejection was simulated in 1 hour time steps along 
the comet’s orbit while it was within 2.5 AU of the Sun. In the second case, ejection 
was simulated to occur at the hour the comet reached perihelion. A 4th order variable 
step-size Runge-Kutta integrator was then used to integrate meteoroid position and 
velocity forward in time, accounting for the effects of radiation pressure, Poynting- 
Robertson drag, and the gravitational forces of the planets, which were computed 
using JPL’s DE406 planetary ephemerides. An impact parameter was computed for 
each particle approaching the Earth to create a flux profile, and the results compared 
to observations of the 1998 and 1999 Leonid showers, and the 1993 and 2004 
Perseids. 
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1. Introduction 

The Marshall Space Flight Center (MSFC) meteoroid stream model simulates 
particle ejection and subsequent evolution from comets in order to provide meteor 
shower forcasts to spacecraft operators for hazard mitigation and mission planning 
purposes. This paper is concerned with simulating the evolution of the Leonid and 
Perseid streams associated with comets 55P/Tempel-Tuttle and 109P/Swift-Tuttle, 
respectively. The model is in a fairly early phase of development, thus the results 
reported here are preliminary. The immediate aim was to compare the peak solar 
longitudes resulting from the model with observations of past Leonid and Perseid 
encounters. 


2. Model 

2.1 Overview 

In modeling particle ejection and subsequent evolution from comets Tempel-Tuttle 
and Swift-Tuttle, the workload is broken into three parts. First the test particles are 
created, then their positions and velocities are integrated forward in time, and finally 
the particles are extracted at specific times of interest. 

At the first step, the JPL HORIZONS Solar System Data and Ephemeris 
Computation Service is used to create cometary ephemerides. Particle state vectors 
are generated for each line in the ephemeris while their physical properties are 
determined from a uniform, random draw on log (3, where (3 is the ratio of radiation 



pressure to the Sun’s gravitation, and an assumed density. Two separate ejection 
schemes (two sets of state vectors) are simulated: Timestep and Perihelion. Ejection 
is simulated in 1 hour time steps along the comet’s orbit while it is within 2.5 AU of 
the Sun, hereafter referred to as the Timestep case. Ejection occurring only at the 
hour of perihelion passage is referred to as the Perihelion case. The model of Jones 
and Brown (1996b) is used to determine particle ejection velocity. The particle 
velocity far from the comet is given by 

K, = 41.7 (sin(«r / 2)) 037 (cos z) 0 ' 519 ,2 m~ u 6 p~ u 3 r h ~ ims 

where a is the semi-angle of the spherical cap of ejection, z is the local solar zenith 
angle, R c is the comet radius, m is the particle mass, p is the particle density, and r h 
is heliocentric distance. Sublimation is taken to occur on the day side of the comet 
within a constrained cap angle (a); for further description of the cap angle and the 
other parameters, the reader is referred to Jones and Brown (1996b) and Jones 
(1995). Other studies (e.g. Brown and Jones, 1998; Gockel and Jehn, 2000; Welch, 
2003) have investigated this ejection velocity model. 

At the second step, a 4th order variable step-size Runge-Kutta integrator is used to 
integrate meteoroid position and velocity forward in time, accounting for the effects 
of radiation pressure, Poynting-Robertson drag, and the gravitational influences of 7 
planets, Venus through Neptune. JPL’s DE406 is used to compute the positions of 
the planets; interpolation is done via cubic spline. Integration of the particle 
ensemble is performed on 9 dedicated AMD Athlon™ 64-bit FX processors. The 
CPU time for each case, as discussed in Section 2.2, is approximately 5 days for the 
Leonids and 2 weeks for the Perseids. 

In the last step, particles are extracted within 0.01 AU of Earth within 1 week 
before and after the expected shower peak. The node-Earth distance for each particle 
is computed, along with an impact parameter (IP) for each meteoroid approaching 
Earth. The impact parameter is defined as follows 

IP = ( r e + h atmos )/ D 

where R E is the radius of the earth, h atmos is the height of the atmosphere, and D is the 
Earth-particle distance at nodal crossing; it is scaled to 1 at the top of the atmosphere. 
It is known that D is not the closest a particle will approach Earth, especially in the 
case of streams associated with low inclination comets. In the case of streams 
associated with high inclination comets, however, the difference between D and the 
actual closest approach distance is inconsequential. Thus the IP is valid for streams 
with large orbital inclinations, as is the case for the Leonids and the Perseids. 

2.2 Inputs 

In modeling the Leonids it was necessary to simulate ejection from comet Tempel- 
Tuttle during the last 12 revolutions, epochs 1600-1965. Three hundred thousand 
particles were ejected during each epoch, in both ejection schemes, for a total of 3.6 
million particles for each case, Timestep and Perihelion. The meteoroid production 
rate was varied with the heliocentric distance as r h ' 5 , based on previous simulations, 
and (3 ranged from ~10' 5 to 10" 2 , resulting in a mass range between approximately 1 
pg and 1 kg, as the density was assumed to be 1000 kg m" 3 . A cap angle of a = 30° 
was chosen. The radius of the comet was assumed to be 2.0 km (Hainaut et al., 
1998). 



For each ejection case and two different cap angles, comet Swift-Tuttle was 
modeled with 900,000 particles ejected over each of its last 9 revolutions dating from 
826 to 1862. This resulted in 4 total Perseid cases with 8.1 million particles each. 
Particle production depended on r h ' 6 (Fomenkova et al., 1995). Again, (3 ranged from 
~10" 5 to 10" 2 , the corresponding masses from 1 kg down to 1 pg, assuming density 
was 1000 kg m’ 3 . Cap angles of 25° and 60° were chosen based on the literature, 
respectively Jones (1995), and Jones and Brown (1996a) and Brown and Jones 
(1998). The radius of the comet was assumed to be 1 1.0 km, based on Boehnhardt et 
al. (1996). 

2.3 Model Comparison 

Per Welch (2003), numerical simulations of meteor showers tend to conform to one 
of two types whereby: 1) large numbers of meteoroids are ejected from a parent 
comet and the subsequent evolution of these particles is followed until the period of 
interest, or 2) the meteoroid orbit with the correct change of period necessary to pass 
through the relevant node at the time of the shower in a given year is determined for 
a given ejection epoch. The MSFC model is of the first type. It is similar to models 
seen in Brown (1999), Brown and Arlt (2000), Gockel and Jehn (2000), and 
Vaubaillon and Colas (2002) in that respect. The models in Kondrat’eva and 
Reznikov (1985), Wu and Williams (1996), Asher (1999), and McNaught and Asher 
(1999) are of the second type. Welch (2003) is a combination of these types. It is 
unknown at this time under which type category the work of Lyytinen (1999) and 
Lyytinen and Van Flandern (2004) falls. 

A brief summary comparing the Leonid MSFC model to the Leonid models of 
other authors is given in Table I. The main categories of comparison involve how 
the cometary ephemerides are created, the ejection velocity law used, the particle 
production dependence on heliocentric distance, the ejection scheme - particles 
ejected along the comet’s orbit or only at perihelion, and the integrator used to 
integrate the position of the comets back in time and/or evolve the particles forward 
in time. The MSFC model uses JPL Horizons to calculate the cometary 
ephemerides, instead of introducing the additional complexity of integrating the 
comets, whose dynamics are significantly influenced by jetting, reflected in the Al 
and A2 terms. Ejection velocity is based on Jones and Brown (1996b) and Jones 
(1995). Type 1 models, with the exception of Vaubaillon and Colas (2002), and 
Welch (2003) have explored this ejection velocity model, among others, with a set 
cap angle of 60°. Vaubaillon and Colas (2002) and Welch (2003) have investigated 
the ejection velocity model of Crifo and Rodionov (1997). The MSFC model’s 
particle production dependence on heliocentric distance is not fixed at r h ' 2 but is 
instead based on observational evidence found in the literature or determined based 
on previous simulations. Two different ejection schemes are compared in this 
model: particle ejection along the comet’s orbit within 2.5 AU and also ejection just 
at the hour of perihelion, like models of Types 1 and 2. In the MSFC model a 4th 
order variable step-size Runge-Kutta integrator is used to evolve the particle orbits, 
like Brown (1999), and Brown and Arlt (2000). Various other numerical integrators 
have been used, including Runge-Kutta-Nystrom 12(10)17 (e.g. Welch, 2003), 
Radau 15 (e.g. Vaubaillon and Colas, 2002), and Stumpff-Weiss (e.g. Gockel and 
Jehn, 2000). 



3. Results and Discussion 

A subset of the total number of particles propagated, those extracted within 0.01 AU 
of Earth within 1 week before and after the expected shower peaks of interest, is now 
examined. As stated previously, an impact parameter (IP) is computed for each 
particle of this set approaching the Earth. In effect, this is the scaled probability that 
the particle will hit Earth. (It must be noted that a stream that approaches Earth, but 
does not appear to cross its path, still has a nonzero probability of striking the Earth.) 

The particle IPs are summed in 0.01° solar longitude bins (this corresponds to a 
time interval of 14.6 minutes), thereby representing cumulative probabilities. 
Ideally, the IPs would be summed in 0.005° bins, corresponding to the Earth moving 
approximately 1 Earth diameter plus the height of the atmosphere. There were not 
enough particles in the subset, however, for this binning scheme, so 0.01° solar 
longitude bins were used. 

A Lorentzian is fit to the binned IP versus solar longitude - essentially the flux 
profile - in order to determine the time of the shower peak. The results are compared 
to observations of the 1998 and 1999 Leonid showers, and the 1993 and 2004 
Perseids. In this way, the goal - to compare the solar longitudes of the main shower 
peak resulting from the model with observations of past Leonid and Perseid 
encounters, thereby validating the MSFC meteoroid stream model’s timing - is 
achieved. 

In the following sections the results of this validation are described. Cross-section 
plots showing the composition of the showers by stream are included and appear as 
parts (a) of Figures 1-4. Each of the particles in parts (a) has an associated IP. These 
IPs are summed in solar longitude bins and, for comparison purposes the amplitudes 
are scaled, thus creating a flux profile, shown in parts (b) of Figures 1-4. The model 
results are simply scaled to the peak ZHR observed so that the time of the observed 
peak and that of the peak predicted by the model can be easily compared. The 
observational data is also shown in parts (b). Finally, predictions from other 
modelers, where possible, are given in tabular form. 

3.1 Leonids 

Figure 1 shows the modeling results of the 1998 Leonids. The shower had two main 
peaks. The first peak, thought to be characterized by old ejecta according to Arlt and 
Brown (1999), was not modeled. The second peak (the nodal peak), however, is 
clearly prominent at solar longitude 235.31 1°+0. 007° (ibid). It is made up of 
particles ejected from recent passages of Tempel-Tuttle, namely 1965, 1932, and 
1899. The MSFC model predicted a peak time at 235.27°+0.01°, a difference of 
about 1 hour. Table II summarizes the results of the 1998 Leonid encounter 
alongside predictions from other authors. 

The 1999 Leonid result can be seen in Figure 2. A storm of activity was observed 
at 235.285°+0.001° (Arlt et ah, 1999). The MSFC model predicted the storm peak, 
consisting mainly of particles from 1899 and 1932 but also including those from 
1965, at 235.282°+0.002°. This is a difference of just under 2 minutes. The MSFC 
results for the 1999 Leonids, and the results of other authors, are presented in Table 
III. 



3.2 Perseids 


Figure 3 illustrates the modeling results of the 1993 Perseids. The main peak 
occurred at 139.53°+0.01° (Bone and Evans, 1996; Rendtel and Brown, 1997). A 
time of 1 39.49 1°±0.002° was determined from the MSFC model, a difference of 57 
minutes. Particles ejected from Swift-Tuttle’s passages in 1862, 1479, and 1610 
were the main contributors. This breakdown is somewhat different from Brown 
(1999). The 1862 stream contributed the most significantly. 

The 2004 Perseid result is shown in Figure 4. The maximum occurred at 
139.443°+0.003° (Arlt, unpublished observations). The MSFC model predicted a 
peak time of 139.42°±0.01°, a difference of 34 minutes. Particles from 1862 made 
up the peak. 

Table IV summarizes the results of the Perseid encounters. 

3.3 Model variables 

Two different ejection cases, Timestep and Perihelion, were simulated for both the 
Leonids and the Perseids. In general, the Perihelion case produces narrower particle 
distributions, as is to be expected. An example of this can be seen for the 2000 
Leonids in Figure 5(a). These two cases yield similar peak times at this preliminary 
analysis stage. It is the Timestep case, however, that is felt to be more physically 
realistic. 

The Leonids were only simulated with one cap angle value. For the Perseids, two 
different cap angles, a=25° and a=60°, were tested. The smaller of the two values 
simulates an exit cone of about 60°; the 60° cap angle corresponds to a hemispherical 
exit cone (Jones and Brown, 1996b). Of the two values, the 60° cap angle model 
generally predicts a solar longitude closer to that of the observed maximum ZHR. 
The 1998 Perseids, shown in Figure 5(b), is a good example of this. With scaling, 
the maximum is clearly distinguishable despite the scatter. The observed peak ZHR 
occurred at solar longitude 139.75°+0.03° (Arlt, 1999). The 60° and 25° cap angle 
models have peaks at 139.72°+0.01° and 139.83°+0.01°, respectively. 

4. Summary 

The MSFC stream model predicts, within about an hour or better, the peak times of 
several Leonid and Perseid encounters. Further refinement is necessary in analyzing 
the results, especially concerning ZHR calculations. Work exploring the effects of 
ejection schemes, the Timestep case and the Perihelion case, and cap angles is 
underway. Additionally, the inclusion of observational data, such as population 
indices, is underway in order to better constrain parameters associated with particle 
ejection from comets. Future work addressing the effects of different cometary 
ephemerides is planned; these effects are expected to be significant. 
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Figure Captions 

Figure 1. Results of the 1998 Leonids, Timestep case, (a) Cross-section in x-y ecliptic coordinates 
showing the composition of the shower by stream. As shown in the legend, the different symbols 
represent various streams by year, and the solid black line represents the Earth’s path. Ejecta from 1899, 
1932, and 1965 played the greatest role in the nodal peak, (b) Comparing the scaled model with 
observations. Observational data was collected by the International Meteor Organization (IMO) and 
published in Arlt and Brown (1999). The model is scaled to the maximum ZHR of the nodal peak. 

Figure 2. Results of the 1999 Leonids, Timestep case, (a) X-Y cross-section plot showing the 
composition of the shower by stream. As shown in the legend, the different symbols represent various 
streams by year, and the solid black line represents the Earth’s path. Revs 2 (1932) and 3 (1899) played 
the greatest role in the storm peak, as well as particles from 1965. (b) Comparing the scaled model with 
observations. Observational data was collected by the IMO and published in Arlt et al. (1999). 

Figure 3. Results of the 1993 Perseids, Timestep case, 60° cap angle, (a) Cross-section in x-y ecliptic 
coordinates showing the composition of the shower by stream. As shown in the legend, the different 
symbols represent various streams by year, and the solid black line represents the Earth’s path. Ejecta 
from years 1862, 1479, and 1610 made up the bulk of the peak, with 1862 being most prominent, (b) 
Comparing the scaled model with observations. Observational data was collected by the British 
Astronomical Association (BAA) and published in Bone and Evans (1996). 

Figure 4. Results of the 2004 Perseids, Timestep case, 60° cap angle, (a) X-Y cross-section plot showing 
the composition of the shower by stream. As shown in the legend, the different symbols represent various 
streams by year, and the solid black line represents the Earth’s path. The maximum was made up from the 



1 862 stream, (b) Comparing the scaled model with observations. Observational data was collected by the 
IMO, and is not yet published (Arlt, unpublished observations). 

Figure 5. (a) X-Y cross-section plot showing the composition of the shower by stream for the 2000 
Leonids. As shown in the legend, the different symbols represent various streams by year, and the solid 
black line represents the Earth’s path. The Perihelion ejection case (in black) is shown superimposed on 
the Timestep case (in gray). Note especially how focused the 1965 and 1932 streams of the Perihelion 
case are, as compared to the Timestep case, (b) Comparing two scaled models with observations for the 
1998 Perseids. The 60° cap angle model (in gray) matches the time of the observed maximum better than 
the 25° model (in black). Observational data was collected by the IMO and published in Arlt (1999). 



Tables 

Table I 


Authors (Leonids) 

Cometary Elements/ 
Ephemerides Source 

Ejection Velocity 

Particle 

Production 

Dependence 

Ejection Scheme 

Integrator 

MSFC 

JPL Horizons 

Jones ejection (Jones, 1995; Jones 
and Brown, 1996b) 

r h ' 5 

Timestep, within 2.5 
AU; Perihelion 

Runge-Kutta 4 

Brown (1999), 
Brown & Arlt (2000) 

Yeomans et al. (1996) 

Jones ejection (Jones, 1995), 
modified Jones and parabolic Jones 
(Brown & Jones, 1998), Crifo (1995) 

rh' 2 

Timestep, within 4 AU 

Runge-Kutta 4 

Gockel & Jehn (2000) 

- 

Jones ejection (Jones, 1995), 
modified Jones and parabolic Jones 
(Brown & Jones, 1998), Crifo (1995) 

r° 

Timestep, within 4 AU 

Stumpff-Weiss 

Vaubaillon & Colas (2002) 

Rocher (2002) 

Crifo & Rodionov (1997) 

r°, weighted 
later by rh’ 2 

Timestep, within 3 AU 

Radau 15 

Kondrat'eva & Reznikov (1985) 

independent 

set at 1 value at perihelion - 0 km/s - 
then iteratively determined 

n/a 

Perihelion 

algorithm of 
Kulikov (1960) 

Wu & Williams (1996) 

n/a 

set at 1 value at perihelion - a mean 
ejection velocity less than 0.6 km/s 

n/a 

Perihelion 

Runge-Kutta-Nystrom 

Asher (1999), 

McNaught & Asher (1999) 

Nakano (1998), 
Yeomans et al. (1996) 

n/a 

n/a 

Perihelion 

Radau 15 

Welch (2003) 

Nakano (1998) 

Jones ejection (Jones, 1995), 
modified Jones (Brown & Jones, 
1998), Crifo & Rodionov (1997) 

Ih' 2 

Timestep, within 3.5 
AU 

Runge- Kutta-Ny strom 






Table II 




McNaught 

Brown 






Year Peak 

& Asher 

& Arlt 

MSFC 

Observ- 

MSFC 




(1999) 

(2000) 

Result 

ations Difference 



1998 Nodal 

- 

- 

235.265 

235.311 

-1.1 hr 



, 1965 

235.26 

235.34 

235.262 

n/a 

n/a 



i 1932 

235.27 

235.22 

235.263 

n/a 

n/a 



t 1899 

- 

- 

235.270 

n/a 

n/a 





Table HI 





Kondrat’eva 

McNaught 






Year 

Peak & Reznikov 

& Asher 

Brown 

Lyytinen 

MSFC 

Observ- 

MSFC 


(1985) 

(1999) 

(1999) 

(1999) 

Result 

ations 

Difference 

1999 

Storm 

See 1899 

See 1899 

See 1899 

235.282 

235.285 

-2 min 

r 

1965 

235.279 

- 

- 

235.272 

n/a 

n/a 

f 

1932 

235.273 

- 

235.270 

235.275 

n/a 

n/a 

l 

1899 235.29 

235.29 

235.3 

235.291 

235.295 

n/a 

n/a 


Table IV 


Year 

Lyytinen & Van 
Flandem (2004) 

MSFC 

Result 

Observations 

MSFC 

Difference 

Main 

Contributors 

1993 

- 

139.491 

139.53 

~57 min 

1862, 1479, 1610 

2004 

139.441 

139.42 

139.443 

-34 min 

1862 


Table I. A brief summary comparing the Leonid MSFC model to the Leonid models of other authors. The 
main categories of comparison involve how the cometary ephemerides are created (most modelers take 
cometary elements from various sources and integrate backwards), the ejection velocity law used, the 
particle production dependence on heliocentric distance, the ejection scheme - particles ejected along the 
comet’s orbit or only at perihelion, and the integrator used to integrate the position of the comets back in 
time and/or evolve the particles forward in time. The top section of the table shows models of Type 1, the 
middle section shows models of Type 2, and the bottom section shows the hybrid. It is unknown at this 
time under which type category the work of Lyytinen (1999) and Lyytinen and Van Flandem (2004) falls; 
these works are not included in this table. 

Table II. This table shows the modeling results of various authors alongside the results of the MSFC 
model and the observations for the 1998 Leonids. The activity peak near the passage of the descending 
node of Tempel-Tuttle is referred to as the nodal peak; the fireball peak made up of older ejecta was not 
modeled. Modeling by McNaught and Asher (1999) and Brown and Arlt (2000) found that the 1965 and 
1932 streams contributed to the nodal peak and calculated separate solar longitudes for each. Binning the 
IPs of ejecta from different epochs separately and fitting a curve to the data, the solar longitudes of the 
maxima for 1965 and 1932 ejecta from the MSFC model are shown for comparison purposes. Material 
from 1899 was also found to be a contributor and shown accordingly. Observations showed a peak at 
solar longitude 235.311 (Arlt and Brown, 1999). An analysis of the subset described in the text, 
combining ejecta from each perihelion passage (1965, 1932, and 1899 included) gives a peak at 235.265, a 
difference of about an hour. Errors are omitted for convenience. 

Table III. This table shows the modeling results of various authors alongside the results of the MSFC 
model and the observations for the 1999 Leonids. A storm of activity was observed at solar longitude 
235.285 (Arlt et al., 1999). In the models of Kondrat’eva and Reznikov (1985), McNaught and Asher 
(1999), Brown (1999), and Lyytinen (1999), as well as the MSFC model, the major (or only) contribution 
to the storm peak came from particles ejected in 1899. Solar longitudes of the maxima for material ejected 
in 1899, 1932, and 1965 are given separately. Binning the IPs of ejecta from different epochs individually 
and fitting a curve to the data, the solar longitudes of the maxima for the 1899-1965 streams from the 
MSFC model are shown for comparison purposes. An analysis of the subset described in the text, 
combining ejecta from each perihelion passage (1965, 1932, and 1899 included) gives a storm peak at 
235.282, a difference of about 2 minutes. Errors are omitted for convenience. 








Table IV. This table shows the modeling results the MSFC model, alongside the results of Lyytinen & 
Van Flandern (2004) where possible, and the observations for the 1993 and 2004 Perseids. The main peak 
of the 1993 Leonids was at 139.53 (Bone and Evans, 1996; Rendtel and Brown, 1997). The MSFC model 
showed a peak at 139.491, a difference of about 57 minutes. Particles ejected from Swift-Tuttle’s 
passages in 1862, 1479, and 1610 were the main contributors. In 2004, material from 1862 was the main 
contributor to the peak. The prediction of Lyytinen and Van Flandern (2004) was very close to what was 
observed. The result of the MSFC model, with a peak at 139.42, was about 34 minutes premature. Errors 
are again omitted for convenience. 




0.55 0.56 0.57 234.0 234.5 235.0 235.5 

x (AU) solar longitude (J2000.0) 


Figure 2 


b) 



0.55 0.56 0.57 235.0 235.2 235.4 235.6 235.8 

x (AU) solor longitude (J2000.0) 


Figure 3 


b) 



0.770 0.775 0.780 139.0 139.2 139.4 139.6 139.8 

x (AU) solor longitude (J2000.0) 





0.55 0.56 0.57 234.0 234.5 235.0 235.5 

x (AU) solar longitude (J2000.0) 


Figure 2 


b) 



0.55 0.56 0.57 235.0 235.2 235.4 235.6 235.8 

x (AU) solor longitude (J2000.0) 


Figure 3 


b) 



0.770 0.775 0.780 139.0 139.2 139.4 139.6 139.8 

x (AU) solor longitude (J2000.0) 








